library(plyr)
library(magrittr)
library(tidyverse)
library(readxl)


rm(list=ls())
home = 'C:/Users/Jason/Dropbox/VNA_Responsiveness/Analysis/JOP-dataverse/'


permutation.results = paste0(home, 'RI-analyses.Rds') %>%
  readRDS

perm.cit = subset(permutation.results, 
                  Result=='tab2col2' & term=='Citizen')$estimate
perm.firm = subset(permutation.results, 
                   Result=='tab2col2' & term=='Firm')$estimate
perm.cf.diff = perm.cit - perm.firm

perm.cit.edu = subset(permutation.results, 
                      Result=='placebo' & term=='Citizen:Education')$estimate
perm.firm.edu = subset(permutation.results, 
                       Result=='placebo' & term=='Firm:Education')$estimate
perm.cf.diff.edu = perm.cit.edu - perm.firm.edu

perm.cit.speak = subset(permutation.results, 
                        Result=='tab3col2' & term=='Citizen')$estimate
perm.firm.speak = subset(permutation.results, 
                         Result=='tab3col2' & term=='Firm')$estimate
perm.cf.diff.speak = perm.cit.speak - perm.firm.speak



experimental.results = paste0(home, 'experimental-analyses.Rds') %>%
  readRDS

exp.cit = subset(experimental.results, 
                 Result=='tab2col2' & term=='Citizen')$estimate
exp.firm = subset(experimental.results, 
                  Result=='tab2col2' & term=='Firm')$estimate
exp.cf.diff = exp.cit - exp.firm

exp.cit.edu = subset(experimental.results, 
                     Result=='placebo' & term=='Citizen:Education')$estimate
exp.firm.edu = subset(experimental.results, 
                      Result=='placebo' & term=='Firm:Education')$estimate
exp.cf.diff.edu = exp.cit.edu - exp.firm.edu

exp.cit.speak = subset(experimental.results, 
                       Result=='tab3col2' & term=='Citizen')$estimate
exp.firm.speak = subset(experimental.results, 
                        Result=='tab3col2' & term=='Firm')$estimate
exp.cf.diff.speak = exp.cit.speak - exp.firm.speak



CI.cf.diff = CI.cf.diff.edu = CI.cf.diff.speak = NULL
for(i in seq(from=-2, to=2, by=0.001)) {
  CI.cf.diff = rbind(CI.cf.diff, 
                     c(shift=i,
                       shifted=exp.cf.diff+i,
                       upper=sum((exp.cf.diff+i)>=perm.cf.diff)/length(perm.cf.diff),
                       lower=sum((exp.cf.diff+i)<=perm.cf.diff)/length(perm.cf.diff)))
  CI.cf.diff.edu = rbind(CI.cf.diff.edu, 
                         c(shift=i,
                           shifted=exp.cf.diff.edu+i,
                           upper=sum((exp.cf.diff.edu+i)>=perm.cf.diff.edu)/length(perm.cf.diff.edu),
                           lower=sum((exp.cf.diff.edu+i)<=perm.cf.diff.edu)/length(perm.cf.diff.edu)))
  CI.cf.diff.speak = rbind(CI.cf.diff.speak, 
                           c(shift=i,
                             shifted=exp.cf.diff.speak+i,
                             upper=sum((exp.cf.diff.speak+i)>=perm.cf.diff.speak)/length(perm.cf.diff.speak),
                             lower=sum((exp.cf.diff.speak+i)<=perm.cf.diff.speak)/length(perm.cf.diff.speak)))
}; rm(i)
CI.cf.diff = as.data.frame(CI.cf.diff)
CI.cf.diff.edu = as.data.frame(CI.cf.diff.edu)
CI.cf.diff.speak = as.data.frame(CI.cf.diff.speak)



CIs.Primary = c((-1)*CI.cf.diff[CI.cf.diff$lower>0.025,][sum(CI.cf.diff$lower>0.025),]$shift, 
                exp.cf.diff,
                (-1)*CI.cf.diff[CI.cf.diff$upper>0.025,][1,]$shift)
CIs.Placebo = c((-1)*CI.cf.diff.edu[CI.cf.diff.edu$lower>0.025,][sum(CI.cf.diff.edu$lower>0.025),]$shift, 
                exp.cf.diff.edu,
                (-1)*CI.cf.diff.edu[CI.cf.diff.edu$upper>0.025,][1,]$shift)
CIs.Speak = c((-1)*CI.cf.diff.speak[CI.cf.diff.speak$lower>0.025,][sum(CI.cf.diff.speak$lower>0.025),]$shift, 
              exp.cf.diff.speak,
              (-1)*CI.cf.diff.speak[CI.cf.diff.speak$upper>0.025,][1,]$shift)



CIs = rbind(CIs.Primary,
            CIs.Placebo,
            CIs.Speak) %>%
  as.data.frame %>%
  set_colnames(c('Lower',
                 'Estimate',
                 'Upper')) %>%
  mutate(Specification=factor(x=c('DV: Preparation\nAll delegates + Covariates',
                                  'DV: Preparation\nAll delegates + Covariates\nAccounting for Hawthorne effect',
                                  'DV: Speaking\nAll delegates +Covariates'), 
                              levels=c('DV: Preparation\nAll delegates + Covariates',
                                       'DV: Preparation\nAll delegates + Covariates\nAccounting for Hawthorne effect',
                                       'DV: Speaking\nAll delegates +Covariates')))
CIs

CIs.plot = ggplot(CIs, 
                  aes(x=Specification, y=Estimate, ymin=Lower, ymax=Upper)) + 
  geom_linerange(size=1) +
  geom_point(size=2) +
  geom_hline(yintercept=0, color='red', linetype=2) +
  scale_y_continuous(breaks=seq(from=-0.3, 
                                to=0.3, 
                                by=0.1), 
                     labels=round(seq(from=-0.3, 
                                      to=0.3, 
                                      by=0.1), 
                                  digits=1)) +
  coord_cartesian(ylim=c(-0.274, 0.274)) +
  labs(x=NULL, y=expression(beta[Citizen]-beta[Firm])) +
  theme_minimal() +
  theme(panel.border=element_rect(color='black', fill=NA))


ggsave(filename='figure-08.png', plot=CIs.plot, path=home, width=9, height=3, units='in')
ggsave(filename='figure-08.tiff', plot=CIs.plot, path=home, width=9, height=3, units='in')
ggsave(filename='figure-08.eps', plot=CIs.plot, path=home, width=9, height=3, units='in', device=cairo_ps)
